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ABSTRACT 

'  Many  multivariate  interpolation  schemea  require  as  data  valuea  of 
derivatives  that  are  not  available  in  a  practical  application,  and  that 
therefore  have  to  be  generated  suitably,  A  specific  approach  to  this  problem 
is  described  that  is  modeled  after  univariate  spline  interpolation. 

Derivative  values  are  defined  by  the  requirement  that  a  certain  functional  be 
minimised  over  a  suitable  space  subject  to  Interpolation  of  given  positional 
data.  In  principle,  the  technique  can  be  applied  in  arbitrarily  many 
variables.  The  theory  is  described  in  general,  and  particular  applications 
are  given  in  one  and  two  variables,  A  major  tool  in  the  implementation  of  the 
technique  is  the  B&zler-Bernstein  form  of  a  multivariate  polynomial.  The 
technique  yields  visually  pleasing  surfaces  and  is  therefore  suitable  for 
design  applications.  It  is  less  suitable  for  the  approximation  of  derivatives 
of  a  given  function. 
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SIGNIFICANCE  AND  EXPLANATION 

This  report  addresses  the  problem  of  the  interpolation  of  bivariate 
scattered  data  ( i)  as  they  arise  for  example  in  the  Computer  Aided 
Geometric  design  of  the  shape  of  a  vehicle  or  in  the  representation  of 
measured  data.  Many  techniques  proceed  by  triangulating  the  domain  and  then 
defining  the  lnterpolant  piecewise  on  each  individual  triangle.  However,  such 
interpolation  schemes  usually  require  as  data  values  of  derivatives  that  the 
user  is  unable  to  supply.  The  report  describes  an  approach  to  the  generation 
of  such  data  from  the  given  information.  It  is  modeled  after  univariate  cubic 
spline  interpolation  and  consists  of  choosing  that  lnterpolant  which  minimises 
a  suitable  functional  (modeling  e.g.  the  strain  energy  of  a  clamped  elastic 
plate)  over  a  suitable  function  space,  subject  to  interpolation.  The  user  has 
soms  choice  in  the  selection  of  the  functional.  The  technique  has  the 
advantage  of  yielding  visually  pleasant  surfaces.  It  has  the  drawback  that 
the  derivative  generation  is  a  global  proceaa  that  depends  upon  all  given 
data.  The  scheme  can  be  considered  either  a  method  of  constructing  a  once 
differentiable  lnterpolant  to  positional  data,  or  a  method  of  constructing 
derivative  values  through  second  order  that  can  then  be  used  to  generate  a 

twice  differentiable  surface  by  employing  another  interpolation  method.  The 

j  Arc'.1 '  1 

basic  ideas  are  also  illustrated  in  the  univariate  context.  !~ 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 
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MULTIVARIATE  SCATTERED  DATA  DERIVATIVE  GENERATION  BY  FUNCTIONAL  MINIMIZATION 

Peter  Alfeld* 

1.  ITROPPCriOW 

Many  multivariate  Interpolation  schemes  require  values  of  derivatives  as 
data  which  are  not  available  in  a  practical  application*  For  examples  of  such 
schemes*  see,  e*g. ,  Barnhill  and  Farin,  1981,  the  papers  by  Alfeld  and  by 
Alfeld  and  Barnhill,  and  other  papers  in  the  proceedings  edited  by  Barnhill 
and  Nielson,  1984,  and  the  finite  element  literature.  Thus,  if  an 
interpolation  problem  is  to  be  solved,  derivative  data  must  be  generated  from 
existing  positional  data. 

In  this  paper,  a  particular  technique  for  derivative  generation  is 
described  which  is  modeled  after  univariate  spline  Interpolation.  In  the 
univariate  case,  it  is  well  known  that  the  interpolating  spline  minimises  a 
variational  principle,  and  it  has  been  recognized  (see  Dube  1975,  Duchon  1976, 
and  Nielson  1983}  that  a  similar  property  is  desirable  also  in  the 
multivariate  context.  Minimising  an  appropriate  functional  over  as  large  a 
space  as  possible  has  the  drawback  that  ususally  this  is  difficult 
computationally  and  that  In  general  the  result  will  not  be  a  function  with  a 
simple  and  convenient  (e.g.  polynomial)  structure.  In  our  approach,  a  simple 
structure  is  forced  by  choosing  the  Interplants  from  a  class  of  piecewise 
polynomial  functions  (defined  over  a  tesselatlon  into  slmplices,  l.e. 

Intervals  or  triangles)  such  that  some  user  specified  functional  is 
minimised.  The  approach  is  related  to  the  finite  element  method  for  the 
solution  of  variational  problems,  and  Indeed  many  computational  techniques 
familiar  in  that  context  can  also  be  used  for  derivative  generation.  The 
technique  is  mathematically  well  founded  and  yields  visually  pleasing 
surfaces*  This  makes  it  suitable  for  design  purposes.  Its  main  drawback  is 

that  the  derivative  generation  Is  a  global  process  that  requires  the 
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eons truct Ion  and  solution  of  a  large  sparse  system  of  linear  equations.  After 
the  preprocessing  stage  of  derivative  generation,  however,  the  schemes  are 
local,  l.e.  the  evaluation  of  the  interpolant  requires  data  only  on  that 
simplex  In  which  the  data  point  resides  . 

The  techniques  presented  here  can  be  viewed  in  at  least  two  different 
ways.  The  point  of  view  that  Is  reflected  In  the  title  Is  that  they  simply 
generate  derivative  values  at  points.  Also  possible  Is  the  view  that  the 
derivative  generation  and  the  interpolation  are  an  entity  that  allows  the  user 
to  construct  smooth  Interpolants  to  positional  data,  without  requiring  user 
specified  derivative  values.  From  that  point  of  view,  the  fact  that  the 
derivative  generation  may  be  divorced  from  the  Interpolation,  and  derivatives 
may  be  used  for  separate  purposes  (e.g.  generating  a  smoother  surface),  is  a 
serendipitous  byproduct  of  the  technique. 

In  traditional  numerical  analysis,  one  usually  thinks  of  an  Interpolation 
scheme  as  a  means  of  approximating  an  underlying  smooth  function.  This, 
however,  Is  often  Inappropriate  in  multivariate  applications.  Rather,  the 
attitude  usually  taken  by  designers  and  engineers  is  that  one  is  given  some 
data  with  the  objective  of  passing  a  visually  pleasing  surface  through  the 
data.  Indeed,  in  a  design  problem,  the  primitive  function  is  the  object  to  be 
created,  rather  than  reproduced  or  approximated.  The  designer  may  even  change 
the  positional  data  In  order  to  get  a  better  surface.  This  Is  why  the  title 
of  this  paper  contains  the  term  "generation"  rather  than  "estimation".  We  are 
Interested  not  in  estimating  derivatives  that  have  no  meaning  in  the  context 
of  the  given  problem,  but  rather  in  trying  to  determine  free  parameters  in  a 
way  that  would  satisfy  a  designer  or  engineer.  This  attitude  Is  the  main 
reaaon  why  the  question  of  errors  is  not  addressed  In  this  paper. 

The  paper  Is  organized  as  follows: 


In  section  2,  vs  devslop  the  theoretics!  frame -work  for  the  technique. 

We  define  the  notation  end  address  the  questions  of  invariance  of  the 
interpolants  under  certain  geometrical  trsns format ions,  of  existence  and 
uniqueness  of  the  minimizing  interpolant,  end  of  the  maximum  degree  of  the 
polynomials  that  are  reproduced  exactly  by  the  method.  This  frame-work 
applies  to  functions  of  arbitrarily  many  variables. 

In  section  3,  ve  consider  the  univariate  case.  We  describe  (piecewise 
cubic  and  qulntlc)  Hermlte  interpolation  and  give  formulas  for  picking  the 
values  of  the  derivatives  so  as  to  minimize  the  norm  of  certain  derivatives  of 
the  Interpolant.  These  include  natural  cubic  and  qulntlc  splines  as  special 
cases. 

In  section  4,  ve  go  on  to  functions  of  two  variables.  This  section  is 
the  heart  of  the  paper.  We  consider  one  basic  Interpolant  (defined  over  a 
triangulation),  namely  a  piecewise  qulntlc  C1  Interpolant  that  Interpolates 
to  gradients  and  Hessians  (as  well  as  to  position).  The  Interpolant  has  been 
long  known  in  the  finite  element  context. 

In  all  cases,  positional  data  are  considered  given,  and  derivative  values 
are  considered  parameters  to  be  picked  so  as  to  minimize  a  functional.  The 
functional  in  all  cases  is  the  Lj  norm  of  some  suitably  picked  set  of 
derivatives. 

In  section  5,  ve  illustrate  the  technique  with  some  bivariate 


r  i* 


2.  THE  F1AHB-W0BK 


la  this  section,  we  define  some  notation  that  will  be  used  In  the 
sequel,  tfe  also  describe  the  approach  In  general,  and  observe  some  simple  but 
fundamental  properties.  The  following  sections  can  then  be  considered 
specific  examples  for  our  general  approach.  Some  of  our  definitions  and 
concepts  could  be  formulated  more  generally,  but  we  restrict  them  to 
situations  as  they  arise  In  practical  applications. 

We  denote  by  D  a  polyhedral  region  in  E**  that  has  been  tesselated 
Into  slmpllces.  We  do  not  address  the  question  of  how  this  tesselation  may  be 
accomplished.  See  Little,  1983,  for  some  answers  to  this  question. 

Our  interpolant  will  be  piecewise  polynomial  (so  that  Integrations  can  be 
carried  out  easily).  We  denote  by  Bm  a  space  of  functions  that  are 

Q 

polynomial  on  each  of  the  slmpllces  in  D  and  that  are  globally  m  times 

continuously  differentiable.  The  space  contains  all  polynomials  of  degree  up 

to  and  Including  n  .  In  one  variable,  B*  will  not  contain  any  higher 

n 

degree  polynomials,  but  in  more  than  one  variable,  B™  may  contain  some 
higher  degree  polynomials,  depending  upon  the  application. 

We  assume  that  our  Interpolant  Is  required  to  interpolate 
to  N  data  defined  by  linear  functionals,  thus,  If  p  denotes  the 

Interpolant , 


(2.1) 


■  d^  ,  1  *  !,«..  ,N 


In  most  applications,  the  will  be  point  evaluations  at  the  vertices 
of  the  tesselation.  We  also  assume  that  the  Interpolant  contains  free 
parameters  p1  that  we  have  to  specify.  In  most  applications,  these 
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parameters  will  be  values  of  derivatives  at  points.  We  express  the  link 
between  the  lnterpolants  and  the  parameters  In  terns  of  linear  functionals : 


(2*2)  P^p  •  •  1  ■  1|... |M  » 


We  collect  the  data  d^  into  the  vector  d  and  now  define  the  efflne 
space  that  contains  ell  candidates  for  the  final  lnterpolant  by 


A*(d)  -  (q  c  B*  :  (2.1)  and  (2.2)  hold  where  p^c  1  } 


Fundamental  to  our  technique  Is  the  minimisation  of  certain 
derivatives.  We  consider  the  differentiation  operator 


ax,  ax.  ...  ax. 


•••»  ^  ®  {lpmmepd)  lj  1 


2 

-)  • 


We  also  use  the  functional  defined  on  Bb  by 

II 


(2.3) 


Fkq  •  if  Dkq 
{T}  T 


Here,  the  summation  Is  over  all  siapllees  In  the  tesselatlon  of  0  and 
the  Integration  Is  over  an  individual  simplex.  For  example.  If  k  •  1  ,  we 
obtain  the  Eudedian  norm  of  the  gradient,  and  if  k  -  2  ,  we  obtain  the 
Frobenlus  norm  of  the  Hessian  of  q  . 

In  the  now  defined  context,  of  the  infinitely  many  lnterpolants  contained 
in  A™(d)  ,  we  choose  that  one  that  minimises  F  *  . 
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The  following  theorems  cover  some  of  the  relevant  properties  of  the  above 
approach. 

U  a 

Theorem  1:  There  exists  a  unique  minlmlser  of  FK  in  A  (d)  for 

n 

all  d  if  and  only  if  the  only  function  in  Aq(0)  piecewise  of 
degree  <  k  is  the  zero  function. 

lammrk:  Thus  we  reduce  the  question  of  existence  and  uniqueness  for  the 
sdnlmizatlon  problem  to  a  similar  question  for  a  much  simpler  Legrange 
interpolation  problem. 

Theorem  2:  If  p  minimizes  F^  over  A®(d)  and  p  is  unique, 

then  p  is  invariant  under  translations,  isotropic  rescalings,  and  rotations. 

Theorem  3:  If  p  minimizes  F^  over  A^(d)  ,  and  p  is 

unique,  k  <  n  ,  and  d^  ■  L^q  for  some  polynomial  q  of  degree  <  k  , 
then  p  ■  q  . 

Proof  of  Theorem  1:  The  functional  F*  is  quadratic  and  minimizing  it 
requires  the  solution  of  a  linear  system.  If  d  ■  0  ,  then  the  linear  system 
la  also  homogeneous  (since  the  zero  function  is  a  solution).  Moreover,  the 
coefficient  matrix  is  Independent  of  d  (as  can  be  seen  by  writing  the 
interpolant  in  cardinal  form).  A  linear  system  with  arbitrary  right  hand  side 
possesses  a  unique  solution  if  and  only  if  the  corresponding  homogeneous 
linear  system  possesses  a  unique  solution.  If  the  zero  function  is  the  only 
function  in  Ar  piecewise  of  degree  <  k  ,  then  the  solution  is  unique  as  any 
Interpolant  of  a  higher  degree  would  yield  a  non-zero  value  for  . 


Proof  of  Theorem  2i  tfo  consider  •  change  of  variables 
T  2 

.  s  ■  Lx+b  where  L4L  •  a  I  for  soae  constant  a  ,  ,  with  I  being  the 

identity  matrix.  This  class  of  transformations  includes  translations 
(L  ■  0),  isotropic  rescaling  (L  ■  al  *  b  ■  0),  and  rotations 
(  a  *  1  .  b  •  0)  .  Ve  have  to  show  that  the  solution  of  the  minimisation 

I 

problem  is  independent  under  this  change  of  variables. 

Let  p(x)  ■  p(Lx+c)  .  We  will  show  thet 

^  (2.4)  D -  a2kDkp  . 


» 


I 

In  the  following,  subscripts  of  p  denote  partial  derivatives,  and  all 

evaluations  take  place  at  Lx+b.  - 

fe 


Thus  the  value  of  the  functional  F^p  is  proportional  to  F^p  which 
establishes  the  theorem. 

To  see  (2.4)  requires  a  simple  but  messy  manipulation  of 

U. 

summations.  D  p  ia  obviously  Independent  of  b  ,  so  we  may  assume  b  ■  0 
We  denote  the  (i,J)  entry  of  L  by  .  We  next  see  by  induction  that 


£ 


k 

IT  L. 


j.  p(Lx  +  b)  ■  I  u  •*.  ,  _  _ 

*i,  a*i,,,*a*i.  jj . Jk  r-1  Vr 


(Lx  +b)  . 


ll 


'1  J2 


| 
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Using  this  relstlon,  we  obtsln: 


Dk;  - 1  .  r  ,  f 

ilti2,...,ik  «  *..ax 


r  l 

uiJi 


..  1*  .  L  .  P 

\Jk  \jk  Hi* 


■  £  .  I  l-  .  .11  4  h  4  .  .  . 

V*\  *1  U1J1  h  “2^2  °2J2 
V*\ 

.  I  L  .  L  ,  p  p  “  <*2k  Ip2  ■  <*2kDkp 

Jk  Vk  V>k  V\  V"\  V  \  Hi*  *  *  *\t. 


Proof  of  Theorem  3:  Suppose  that  the  data  are  generated  by  a 

B  lr 

polynomial  q  of  degree  <  k  <  n  .  Then  q  e  AQ(d)  and  q  assumes  the 
minimum  possible  value  zero.  Since  the  mlnimlzer  Is  unique  we  have  that 
P  ■  q  • 

In  the  sequel,  we  will  make  use  of  a  technique  familiar  from  the  finite 
element  technique  (see,  e.g.,  Strang  and  Fix,  1973,  p.28).  Rather  than 
integrating  over  the  entire  domain  D  and  then  differentiating  with  respect 
to  the  parameters,  we  Integrate  over  an  individual  simplex  (element), 
differentiate  with  respect  to  those  parameters  that  contribute  nontrlvlally  on 
the  element,  and  In  this  fashion  obtain  an  element  stiffness  system.  The 
element  stiffness  systems  then  have  to  be  assembled  into  the  global  stiffness 
system,  and  that  system  has  to  be  solved.  We  assume  the  assembly  technique  to 
be  familiar.  Also,  we  will  not  address  the  question  of  how  the  (large, 
sparse,  positive  definite)  global  stiffness  system  may  be  solved  except  to 
point  out  that  this  is  usually  done  by  direct  methods  and  that  sophisticated 
software  for  this  purpose  exists  (see  Ceorge  and  Liu,  1981). 


3.  WIVABIATE  APPLICATIONS 


3.1  Th«  C*  Cage 


Suppose  we  ere  given  decs  (xn,fn)  ,  n  -  ,  where 


•  -  xx  <  *2  <  ...  <  xN  -  b 


We  denote  the  Interval  by  write  ^  •  xn+l"xn  (*or 
n  ■  1,2,... ,N-1)  .  Our  objective  is  to  interpolate  to  the  given  data,  l.e.  we 
are  looking  for  a  function  s  ,  say,  that  satisfies 


(3.1) 


a(xn)  ■  fn,  for  all  n  ■  1,2,...,N 


In  this  subsection,  we  require  s  to  be  in  the  space  B^  ,  ,  l.e.  the 
space  of  all  functions  that  can  be  represented  as  a  cubic  polynomial  on  each 
subinterval  IQ  ,  and  that  are  once  continuously  differentiable  on 
I:  “  (a,b].  In  the  notation  of  section  2,  we  let  d^-f^  and 
Pj  ■  p'(x1)  where  p  is  the  interpolant.  The  affine  space  is  the 

space  of  all  functions  in  ,  that  satisfies  the  interpolation  condition 
(3.1).  The  functional  (for  k  -  0,1, 2, 3)  becomes 


(3.2) 


F-P  -  V  L1+1  (p<»)2 


Jx 

1-1  i 
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l 


To  sake  the  notation  of  the  following  reaulta  more  compact,  we  consider  the 
following  more  general  functional 


k  k 

G  P  -  I  w  F  P 
k-0 


where  the  are  weight  factors  that  can  be  chosen  by  the  user.  For 

Instance,  the  choice  w2  "  1  •  w0  "  W1  "  w3  ”  0  yields  the  natural 

lnterpolatory  cubic  spline  (which,  remarkably,  is  twice  continuously 
differentiable  everywhere  on  I  ). 

The  differentiation  and  integration  on  an  individual  subinterval  was 
carried  out  using  the  symbol  manipulation  language  REDUCE  (Hearn,  1983).  We 
obtain  the  element  stiffness  matrix 


A  -  c 
n 


*11  *12 
*21  *22 


and  the  right  hand  side  gQ  ■  c 


*  J 


where 


210  h 


5 


l  J 
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(3.3) 


•ll  ■  *22  ■  *hn"0  +  56hn,'l  +  ,680!'^"2  +  15120  *3 
*12  •  *21  -  ”3hnW0  '  1‘howl  +  **°^nv2  +  15120  *3 


rl 


r2 


■  £  i-"‘<22Vi  +  13I„>»0  +  ‘2»‘<f 
+  +  502*°<  Vl 

*  £  +  1Sf.>»0  +  421«i<f 

♦  2520hn(Vr'a>»2  +  30240“»H 


n+l*fn)y 


-f.)y31 


n+1 


-f  )w 
n 


-fn)y3> 


1 


The  formulas  given  in  (3.2)  and  (3.3)  can  be  assembled  into  a  global 
stiffness  ayatem  vhoae  aolutlon  deflnea  a  piecewise  cubic  Hermlte  lnterpolant 
that  minimises  the  weighted  derivatlvea  defined  in  (3.2). 

Of  apeelal  interest  is  the  case  where  w^  ■  1  and  the  other  weights 
equal  aero.  The  element' stiffness  matrix  is  of  rank  one.  The  assembly 
procedure  yields  a  consistent  system  of  equetlons  with  rank  deficiency  one. 
Thus  there  is  a  one  parameter  family  of  minimizing  interpolants.  That  this  is 
so  can  also  be  recognized  directly:  Any  piecewise  quadratic  lnterpolant  would 
minimize  the  norm  of  the  third  derivative.  A  quadratic  polynomial  has  three 
degrees  of  freedom.  On  the  first  subinterval,  we  have  to  Interpolate  to  two 
function  velues,  retaining  one  free  parameter.  On  each  subsequent  interval, 
we  have  to  Interpolate  to  two  function  values,  and  also  have  to  match  the 
derivative  of  the  quadratic  on  the  left  neighboring  Interval,  using  up  all 


» 


9 


9 


I 
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three  parameters.  Thus,  overall  one  parameter  remains  at  our  disposal.  The 
significance  of  this  observation  is  that  quadratic  precision  of  a  cubic 
Hermite  scheme  cannot  be  achieved  by  specifying  the  functional  (3.2)  alone. 
Rather,  an  additional  condition  has  to  be  Imposed  to  eliminate  the  free 
parameter.  A  simple  approach  would  consist  of  determining  the  first 
derivative  at  some  grid  point  by  using  a  quadratically  precise  interpolation 
formula,  and  then  requiring  the  minimizing  interpolant  to  reproduce  that 
derivative. 


3.2  The  C2  Came 

The  development  in  this  subsection  parallels  that  in  the  preceding  one, 

2 

except  that  we  consider  interpolants  in  ,  and  also  consider  second 

derivatives  of  the  interpolant  at  the  grldpoints  as  parameters.  The 

parameters  at  our  disposal  are  dQ  -  p' (xR)  and  dQ|i;  -  p"(xn)  for 

2 

n  ■  1,...,N  .  The  affine  space  is  the  space  of  all  functions 
2 

in  Bj  that  satisfy  the  interpolation  condition  (3.1). 

Similarly  as  before,  we  consider  the  general  functional 

JS>  "  E 

k-0  * 

L 

where  the  F  are  defined  in  (3.2).  Carrying  out  the  integration  and 
differentiation  on  a  subinterval  ZQ  and  writing  h  ■  hjj  , 

df  -  f  ..  -  f  ,  p'  -  p'(x  )  ,  p"  -  p“(x  ) 
nTi  n  n  n  n  n 

yields  the  element  stiffness  equation: 
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The  following  cases  ere  of  special  interest: 

w^  ■  6^ ;  The  element  stiffness  matrix  is  non-singular.  The  unique 

solution  for  H  •  1  is  the  linear  lnterpolant.  In  general,  the  interpolant 

2 

minimises  the  norm  of  the  second  derivative  over  the  space  which  contains 

the  interpolating  natural  cubic  spline.  Since  the  spline  minimises  the  norm 

of  the  second  derivative  over  the  larger  space  containing  all 
2 

piecewise  C  functions  that  are  once  differentiable  at  the  grid  points,  the 
piecewise  qulntic  Hermits  interpolant  generated  here  actually  is  the 
Interpolating  natural  cubic  spline. 

w^  -  4^:  The  element  stiffness  matrix  is  of  rank  3  .  In  the  case  of 

N  ■  2  any  Interpolating  quadratic  will  minimise  the  functional  .  In  the 

case  N  ■  3  ,  the  Interpolant  is  the  unique  interpolating  quadratic  which 
2 

renders  C  sero.  The  global  stiffness  matrix  is  non-singular  whenever 
N  >  2  .  In  that  case,  the  Interpolant  is  the  natural  qulntic  spline,  which, 
remarkably,  is  four  times  continuously  differentiable  everywhere  in  I  . 

w^  “  6^:  The  global  stiffness  matrix  is  rank  deficiency  2  for  any  value 
H  >  1  .  Any  piecewise  cubic  Interpolant  will  render  sero. 

w^  ■  6^:  The  global  stiffness  matrix  is  of  rank  deficiency  W-3  for 
any  N  >  1  .  Any  piecewise  quartlc  interpolant  will  render  the  value  of  H*1 


B IV  AIT ATE  APPLICATIONS 


4. 1  Barycentric  Coordinates  «ad  the  Sexier  form  of  a  Polynomial 

Conceptually,  we  proceed  as  in  the  unlvarlaee  case*  However,  the 
algebraic  aanipulations  grow  much  sore  complicated.  The  problem  become* 
tractable  only  by  uelng  barycentric  coordinate*  and  the  Besier  form  of  e 
polynomial.  Ue  develop  these  concepts  for  the  bivariate  case.  The 
generalisation  to  more  than  two  variables  is  straightforward. 

He  restrict  our  attention  to  a  general  triangle  T  with  vertices 
Vj  ,  Vj  ,  and  V3  ,  and  express  the  location  of  a  point  P  in  T  by 

3  3 

(4.1)  P  •  I  b.V.  where  £  b.  ■  1 

i-1  11  i-1  x 

The  are  the  barycentric  coordinates  of  P  .  They  are  determined 
uniquely  by  the  linear  system  (4.1)  provided  the  underlying  triangle  T  is 
non-degenerate  (which  we  will  assume  throughout). 

We  write  a  polynomial  p  of  degree  n  in  the  Besier  form 


(4.2) 


r+«+t-n 


nl  .r.s.t 

rlsit!  crst  blb2b3 


The  computation  of  the  element  stiffness  system  requires  the  following 
two  Ingredients: 

* 

1.  Spatlsl  differentiation  of  the  Bezier  form  of  a  polynomial. 

2.  Integration  of  the  Bezier  form  over  a  general  triangle. 


Ue  address  each  ingredient  in  turn: 


4.1,1  Spatial  Differentiation  of  the  Bezier  form  of  a  Polynomial. 

a 

Ue  let  0  denote  a  spatial  differentiation  operator,  l.e.  V  m  —  for 
some  direction  e  .  Ue  consider  the  general  Bezier  form  given  by  (4.2). 
Then  Pp  is  a  polynomial  of  degree  n-1  which  can  be  written  in  the 

form 


Pp 


I 

r+a+t*n-l 


(n-l)f  - 
rlslt!  rst 


It  can  be  easily  verified  (see  Farln,  1980)  that 


In  our  Application,  wo  oro  interested  in  differentiating  p  several 
tlaas  in  different  directions.  For  later  reference  we  list  the  relevant 
formulas.  Let  V  ,  1  ■  ,  denote  spatial  directional  derivatives 

(such  as  the  partial  derivatives  with  respect  to  x  and  y).  Then 


Vk-1-V " 


r+s+t-n-k 


(n-k)i 

rislt! 


c<k>  b,b,b, 
rst  1  2°3 


Micro 


(4.3) 


and 


c 


rst 


(4.4) c 


rst  )  "  (o“1)lcr+i,s,tPi+lbl  +  cr,s+l,tPi+lb2  +  cr,s,t+iri+l  b3J 


for  1  m  0,1,... ,k  . 

4.1.2  Integration  of  the  Bezier  form  over  a  General  Triangle 

Let  p  be  defined  by  (4.2).  The  polynomial  is  defined  in  terns  of 


barycentrlc  coordinates,  but  integration  is  required  in  terns  of  cartesian 
coordinates.  A  simple  change  of  variables  (see,  e.g.,  Lang,  1968,  p.421) 
yields  for  any  lntegrable  function  f 


4.2  4 


belli  lnterpolent 


The  underlying  lnterpolent  here  le  the  Besler  fora  of  Qtg  ,  e  reduced 
quintlc  lnterpolent  requiring  function  value* ,  gradients,  and  Hessians  at  the 
vertleea  of  I  .  The  overall  piecewise  quintlc  lnterpolent  is  continuously 
differentiable  everywhere  in  the  triangulation.  Tor  the  derivation  of  an 
explicit  expression  for  Q^g  see  Barnhill  and  Ferin,  1981.  For  eoapleteness 
and  reference  we  list  the  relevant  formulas ,  also  introducing  the  notation 
used  in  the  sequel.  The  lnterpolant  is  given  by: 


(4.5) 


a  .  51  . 

Q  rjt*  rl*"'  *•*  121 


where: 


(4.6) 


d500  “  Pl»  d050  “  P2  *  d005  *  P3 


d410  *  5  P3,l  +  d500  *  d401  “  5  P2,l  *  d500  *  d140  “  5  P3,2  *  d050 


(4.7) 


d041  *  5  Pl,2  *  d050  *  d014  *  5  Plf3  *  d005  •  d104  *  5  P2,3  +  d005 


d320  "  20  P33,l  +  2d4l0  “  d500  ’  d230  +  20  P33,2  +  2d140  “  ^50 


d032  *  20  P11 ,2  +  2d041  ”  d050  *  d023  “  20  Pll,3  +  2d  014  “  d005 
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+  »l.a  -  <!,3  '  3),33.2  *  2,13.2' 

• 

i 

and: 

* 

• 

(4.11) 

1 

- 

Vt  -  Q(V1) 

_ 

1 

(4.12) 

Vi  -  <v 

• 

(4.13) 

i 

f  m  — L2 (V  ) 

rJk,i  e^ 

9 

(4.14) 

T 

i 

*i+l*i 

9 

1 

*i  T 

Vi 

■  | 

t  ’ 

and  e  •  V1+1-Vt 

(where  index  arithmetic  is  modulo  3). 

Thus,  the  data  In  (4.11)  are  given  positional  values,  the  parameters  in 

► 

(4.12)  and  (4.13)  are  derivatives  (of  the  global  lnterpolant)  at  vertices 

which  are  to  be 

determined  so  as  to  minimize  a  suitable  functional  defined 

below,  and  the  parameters  defined  in  (4.14)  are  projections  of  edges  onto 

i 

other  edges  that 

ensure  global  differentiability  of  the  lnterpolant. 

• 

It  is  Important  to  keep  in  mind  that  ultimately  we  wish  to  generate 

derivative  data 

(i.e.  gradients  and  Hessians)  at  all  vertices  of  a  given 

i 

trlangulstlon. 

The  local  lnterpolant  Q  is  most  conveniently  expressed  in 

• 

— 

i 
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* 

*  J 

terms  of  derivatives  in  the  direction  of  edges  (which  would  be  inappropriate 
on  a  neighboring  triangle  utilising  common  vertex  data).  Hence  we  need  to 
express  edge  derivatives  in  terms  of  partial  derivatives  in  x  and  y  .  This 
la  accomplished  by  writing: 

(4.15)  gi  -  7)(V1),  -  72Q(Vi) 


and 


T  T 

(4.16)  P  “eg,  P  ■  e  H  e 

k,i  k  i  jk.i  J  i  k 

We  now  proceed  as  in  section  3.  Ve  consider  the  space  of  all 
functions  that  can  be  expressed  as  Qjg  on  each  triangle,  and  that  are 
globally  once  continuously  differentiable.  (Notice  that  this  set  contains  all 
quartlc  polynomials,  but  only  some  quint ic  polynomials).  The  affine 
space  of  Interest  is  the  set  of  all  functions  s  in  that  satisfy  the 
interpolation  requirement 


(4.17)  s(x  ,y  )  “  f  (»  P  ,  j  e  {1 ,2 ,3 }  on  a  particular  triangle) 

i  i  i  J 

We  want  to  minimize  the  functional  defined  in  (2.3)  over  .  The 
parameters  at  our  disposal  are  the  values  of  gradients  and  Hessians  at  the 
grldpoints  (xj^.yj)  . 
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Zn  addition  to  tha  item*  covered  in  ganaral  in  taction  4.1,  va  naad  to 
diffarantiate  tha  lntarpolant  with  respect  to  tha  parameters.  At  bafora,  wa 
conaider  a  tingle  element ,  i.e.  a  tingle  triangle.  On  each  triangle,  va  have 
15  parameters  at  our  ditpoaal.  Va  denote  by  v  ■  [v^]  tha  vector  of  theta 
parameter!.  Bacauaa  of  the  tquara  under  tha  integral  in  (2.3),  tha 
contribution  *(v)  ,  aay  to  F*  from  a  tingle  triangle  la  a  quadratic 
function  of  v  .  Thia  implies  that  we  can  wrlta 


1  T  T 

$(v)  •  —  vAv  +  gv  +  c 

2 

where  A  is  a  symmetric  positive  semidefinite  15  by  15  matrix  with  entries 
•jj  ,  g  it  a  vector  with  entries  ,  and  c  it  a  constant.  The 
contribution  of  $(v)  to  the  global  stiffness  system  is  determined  by  A  and 
g  .  Since  A  it  the  Hessian  of  $(v)  we  can  write 


(4.19) 


i  ,  A  ■ 

ij 


Once  the  a^  have  been  determined,  the  g  can  be  computed  by 
substituting  the  unit  vectors  for  v  ,  thus 


m 


9 


9 


9 


9 
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(where  e^  is  the  unit  vector  with  1  in  the  i-th  place  and  zeros 
elsewhere,  and  c  ■  4(0)).  It  is  therefore  crucial  to  be  able  to  carry  out 
the  differentiation  (4.19).  Applying. the  product  rule  twice  yields 


The  differentiation  of  Q  with  respect  to  the  entries  of  v  can  be  carried 
out  by  differentiating  in  the  recursion  formulas  (4.3)  and  (4.4). 
Differentiating  in  (4.4)  yields 


3  d(i+l) 
rst 


<5-‘> 


+  _LdU> 

3v^  r,s+l»t 


D 


i+lb2 


+  d(1> 

9va  r,s,t+l 


°i+lb3J 


The  recursion  is  initialized  by  differentiating  in  (4.3).  The 
derivatives  can  be  derived  easily  from  the  recursions  (4.6)  through  (4.16). 
However,  the  situation  is  complicated  by  the  fact  that  the  lnterpolant  is 
expressed  in  terms  of  barycentric  coordinates,  whereas  the  derivatives 
entering  as  data  are  cartesian  derivatives.  Computational  efficiency  can  be 
gained  by  observing  that  most  of  the  derivatives  obtained  from  (4.4)  are  zero, 
and  by  avoiding  having  the  program  multiply  by  the  zero  derivatives. 
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For  convenience,  the  non-sero  derivatives  ere  listed  in  the  appendix. 
They  vere  computed  using  the  symbol  manipulation  language  REDUCE  (see  Hearn, 
1973).  In  order  to  avoid  typing  and  typs'settlng  errors,  the  REDUCE  output 
has  been  reproduced  photographically. 


5.  BIVARIATE  EXAMPLES 


In  this  section,  we  illustrate  the  technique  described  in  the  preceding 
sections  by  a  simple  bivariate  example*  In  a  typical  application,  no 
underlying  primitive  function  will  be  available.  However,  for  the  sake  of 
illustration  we  use  a  primitive  first  introduced  by  Franke,  1982.  The  surface 
exhibits  several  phenomena  (two  maxima,  one  minimum,  and  one  saddle  point) 
that  may  occur  in  surface  defined  by  a  true  data  set.  The  primitive  function 
is  given  by 


.  S  .-« W  +  <9j-2>2>/* 

.3  -(<9x+1>2/49  +  (9y+l)/10) 

*  4  * 

.  1  ~((9x-7)2  +  <9y-3)2)/4 

7  * 

_  I  -((9x-4)2  +  (9y-7)2) 

5  * 


Approximations  of  f  on  the  unit  square  were  constructed  from  positional 
values  at  only  36  points.  Figure  1  shows  a  contour  plot  of  f  together  with 
the  data  points  and  the  triangulation  of  the  data.  The  triangulation  was 
constructed  using  Little's  method  described  in  Barnhill,  1977. 

Figure  2  shows  a  hidden  line  plot  of  f  .  All  subsequent  plots  are 
drawn  as  viewed  from  the  same  point  and  with  the  same  vertical  scale.  Figures 
3  through  6  show  the  interpolants  obtained  by  minimizing  the  first  through 
fourth  derivatives  respectively.  Minimizing  second  or  third  derivatives 


yields  a  very  pleasing  surface.  The  interpolant  is  somewhat  shallower  than 
the  primitive  function.  This  is  due  to  the  fact  that  none  of  the  data  points 
is  placed  directly  at  the  extrema,  (cf.  Figure  1.)  which  are  therefore 
poorly  represented  by  the  lnterpolants. 

Minimising  the  gradient  yields  a  surface  that  Is  unlikely  to  be  useful 
for  design  applications.  However,  it  remarkably  resembles  a  scene  one  might 
encounter  in  a  physical  mountain  range.  It  is  Intriguing  to  speculate  that 
minimizing  the  gradient  may  sometimes  be  appropriate  for  generating 
topographical  maps  since  erosion  effects  are  the  harsher  the  steeper  the 
terrain.  Minimizing  fourth  derivatives  generates  large  undulations  and  is  not 
a  good  general  purpose  method. 
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Primitive  Function  -  Contours  and  Triangulation 


FIGURE  1 


Primitive  Function 


ttiim 


m 


m 


FIGURE  2 


Minimization  of  2nd  derivatives 


Minimization  of  3rd  derivatives 


Wm 
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Minimization  of  4th  derivatives 


FIGURE  6 


OOBCLDSIOWS 


A  technique  for  derivative  generation  hat  bee  proposed.  It  is 
particularly  useful  in  a  multivariate  context.  Its  main  properties  are: 

1.  The  preprocessing  stage  of  derivative  generation  requires  the  construction 
and  solution  of  a  large  linear  systea.  However,  this  is  closely  related  to 
the  finite  element  method,  and  existing  well  developed  computational 
techniques  can  be  employed. 

2.  Once  derivatives  have  been  generated,  the  lnterpolant  is  local. 

3.  The  formulas  presented  here  apply  to  the  univariate  and  the  bivariate  ease 
only.  However,  conceptually  the  extension  of  the  technique  to  other  basic 
lntcrpolants  is  straightforward.  The  particular  interpolants  considered  in 
this  paper  are  piecewise  polynomial  functions.  This  Is  convenient  since  the 
technique  calls  for  the  integration  of  functions  with  a  similar  structure  as 
the  lnterpolant.  In  principle,  interpolants  with  a  non-polynomial  structure 
can  also  be  used.  They  would  necessitate  either  the  development  of  special 
integration  formulas  or  the  use  of  numerical  (multivariate)  integration. 

A.  Primarily,  our  technique  is  intended  to  be  used  as  a  means  of  constructing 
visually  pleasing  surfaces.  The  special  application  described  in  section  3. 
can  be  considered  a  bivariate  C*  interpolation  scheme  that  requires  only 
C°  data.  This  is  remarkable  in  view  of  the  fact  that  most  existing  local 
multivariate  interpolation  scheme  require  data  of  at  least  the  same  degree  of 
ssBothnese  as  the  degree  of  smoothness  of  the  lnterpolant.  However,  the 
scheme  also  supplies  values  of  derivatives  at  certain  points.  These 
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derivative  value*  can  be  used  for  other  purposes.  For  example,  the  gradients 
and  Hessians  generated  by  minimising  our  plapewlse  quintlc  C*  lnterpolant 
can  be  used  as  data  for  the  second  of  the  two  piecewise  rational  C* 
lnterpolants  described  In  Alfeld,  1984a.  That  lnterpolant  would  of  course  not 
possess  the  minimisation  property  of  the  polynomial  lnterpolant. 
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out  using  the  symbol  manipulation  language  REDUCE  (Hearn,  1983).  The  Figures 
were  generated  using  the  graphics  packages  PLOT79  (Beebe,  1980).  The  College 
of  Science  of  the  University  of  Utah  provided  computer  time. 
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APPENDIX 


Derivatives  of  the  coefficients  of  Q  with  respect  to  the 
derivative  parameters 

The  following  list  is  computer  generated  which  accounts  for  the  awkward 
notation: 


T  2  Qxx  Qxy 

Q  “  lQ*.Qyl  ,  Q  -  , 

Qxy  Qyy 

a  •  si,  s  -  s2,  a  -  s3 

1  2  3 

T  T  T 

e  •  (elx.ely]  ,  e  ■  [e2x,e2y]  ,  e  ■  [e3x,e3y] 

1  2  3 


d  -  d211,  d  -  dl21 ,  d  -  dll2 
211  121  112 


Only  non-zero  derivatives  are  listed. 


differentiation  with  respect  to  Qx(vl): 
d212:  (5*e2x*s2  -  4*e2x  +  e3x)/10 

d221 :  (  -  e2x  +  5*e3x*s3  -  e3x>/10 

d302:  (  -  2*e2x)/5 

<1311 :  (  -  e2x  +  e3x)/5 

d320:  <2*e3x)/5 

d401:  (  -  e2x)/S 

e3x/5 


d410 


differentiation  with  respect  to  Qy(vl): 
d2l2:  (5*e2y*s2  -  4*e2y  +  e3y)/10 

d221s  (  -  e2y  +  5*e3y*s3  -  e3y)/10 

d302;  (  -  2*e2y)/5 

d3ll:  (  -  e2y  +  e3y)/5 

d320:  (2*e3y)/5 

d401:  (  -  e2y)/5 

d410:  e3y/5 

differentiation  with  respect  to  Qxx(vl): 

d2l2:  (e2x*(  -  5*e2x*s2  +  3*e2x  -  2*e3x))/60 

d221:  (e3x*(  -  2*e2x  +  5*e3x*s3  -  2*e3x))/60 

d302:  e2x**2/20 

d3 11:  (  -  e2x*e3x)/20 

d320:  e3x**2/20 

differentiation  with  respect  to  Qxy(vl): 

d212s  (  -  5*e2x*e2y*s2  +  3*e2x*e2y  -  e2x*e3y  -  e2y*e3x)/30 

d221:  (  -  e2x*e3y  -  e2y*e3x  +  5*e3x*e3y*s3  -  2*e3x*e3y)/30 

d302:  (e2x*e2y)/10 

d311:  (  -  (e2x*e3y  +  e2y*e3x))/20 

d320:  (e3x*e3y)/10 


differentiation  vlth  reepect  to  Qyy(vl): 

d212s  (e2y*(  -  5*e2y*a2  +  3*e2y  -  2*a3y))/60 

d221s  (e3y*(  -  2*a2y  +  5*e3y*s3  -  2*e3y))/60 

d302:  e2y**2/20 

d3U:  (  -  e2y*e3y)/20 

d320:  e3y**2/20 

differentiation  with  reepect  to  Qx(v2): 
d032:  (2*elx)/5 

d041:  elx/5 

dl22r  (5*elx*al  -  elx  -  e3x)/10 

dl31s  (elx  -  e3x)/5 

dl40:  (  -  e3x)/5 

d22ls  (elx  +  5*e3x*s3  -  4*«3x)/10 

d230:  (  -  2*e3x)/5 

differentiation  with  respect  to  Qy(v2): 
d032s  (2*eiy)/5 

d041 i  ely/5 

dl22:  (5viely*sl  -  ely  -  e3y)/10 

dl31:  (ely  -  e3y)/5 

dl40s  (  -  e3y)/5 

d221s  (ely  +  5*e3y*a3  -  4*e3y)/10 

d230:  (  -  2*e3y)/5 
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differentiation  with  respect  to  Qxx(v2): 
d032:  elx**2/20 

dl22:  (elx*(5*elx*sl  -  2*elx  -  2*e3x))/60 

dl31:  (  -  elx*e3x)/20 

d221:  (e3x*(  -  2*elx  -  5*e3x*s3  +  3*e3x))/60 

d230:  e3x**2/20 

differentiation  with  respect  to  Qxy(v2): 
d032:  (elx*ely)/10 

dl22:  (5*elx*ely*sl  *  2*elx*ely  -  elx*e3y  -  ely*e3x)/30 

dl31:  (  -  (elx*e3y  +  ely*e3x))/20 

d221:  (  -  elx*e3y  -  ely*e3x  -  5*e3x*e3y*s3  +  3*e3x*e3y)/30 

d230:  (e3x*e3y)/10 

differentiation  with  respect  to  Qyy(v2): 
d032:  ely**2/20 

dl22:  (ely*C5*ely*sl  -  2*ely  -  2*e3y))/60 

dl31t  (  -  ely*e3y)/20 

d221:  (e3y*(  -  2*ely  -  5*e3y*s3  +  3*e3y))/60 

d230:  e3y**2/20 
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differentiation  with  reepeet  to  Qx(v3): 

d014:  (  -  elx)/5 

d023:  (  -  2*elx)/5 

dl04:  e2x/5 

dll3:  <  -  elx  +  e2x)/5 

dl22:  (5*elx*el  -  4*elx  +  e2x)/10 

d203s  (2*e2x)/5 

d212:  (  -  elx  +  5*e2x*e2  -  e2x)/10 

differentiation  with  reepeet  to  Qy(v3): 

d014i  (  -  ely)/5 

d023:  (  -  2*ely)/5 

dl04:  e2y/5 

dll3:  (  -  ely  +  e2y)/5 

dl22:  (5*ely*el  -  4*ely  +  e2y)/10 

d203:  (2*e2y)/5 

d212s  (  -  ely  +  5*e2y*e2  -  e2y)/10 

differentiation  with  reepeet  to  Qxx(v3): 
d023:  elx**2/20 

dl 13s  (  -  elx*e2x)/20 

dl22s  (elx*(  -  5*elx*el  +  3*elx  -  2*e2x))/60 

d203:  e2x**2/20 

d212:  (e2x*(  -  2*elx  +  5*e2x*s2  -  2*e2x))/60 


differentiation  with  respect  to  Qxy(v3): 

4023:  <elx*ely)/10 

dl 13 :  (  -  (elx*e2y  +  ely*e2x))/20 

dl22:  (  -  5*elx*ely*sl  +  3*elx*ely  -  elx*e2y  -  ely*e2x)/30 

d203:  (e2x*e2y)/10 

4212:  (  -  elx*e2y  -  ely*e2x  +  5*e2x*e2y*e2  -  2*e2x*e2y)/30 

differentiation  with  reapect  to  Qyy(v3): 

4023:  ely**2/20 

dl!3:  (  -  ely*e2y)/20 

4122:  (ely*(  -  5*ely*el  +  3*ely  -  2*e2y))/60 

d203:  «2y**2/20 

4212:  (e2y*(  -  2*ely  +  5*e2y*a2  -  2*e2y))/60 
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